House Price - Modelagem

1 Objetivos

  • Quais Modelos utilizar ?
  • Modelos no Radar
    • Regressão por Redes Bayesians
    • Modelos Lineares Generalizados
    • Modelos de Machine Learning
  • Métricas de Performance
    • RMSE

2 Conjunto de dados

model.train <- fread('../outputs/model.train.csv', 
                     sep=",", 
                     showProgress = FALSE)[,-1] %>%
                     data.frame(stringsAsFactors = T) %>% 
                     select(Id,SalePrice,everything())  
tipo <- lapply(model.train,class)
model.train[,unlist(tipo) != 'integer'] <- data.frame(apply(model.train[,unlist(tipo)!='integer'],2,factor))
model.train

2.1 Importancia das variáveis

Obs: Esta parte do código ficará comentada para o markdown rodar mais rápido e para ilustrar o resultado irei apenas colocar a imagem da sugestão de seleção de variáveis feita pelo boruta, caso queiram rodar o código é só descomenta-lo.

selecao-de-variaveis

selecao-de-variaveis

3 Ajustando um modelo aos dados

Como a distribuição da variável SalePrice é assimétrica com cauda a direita, podemos utilizar um MLG, modelos linerares generalizados, em particular a distribuição Gamma.

Em um primeiro momento utilizei as variáveis selecionadas pelo Boruta, 50 variáveis, porém após essa pré-seleção observei que algumas delas estavam com os p-valores muito alto então apliquei o método StepAic para continuar o processo de seleção o que me levou a configuração final abaixo.

Ainda é importante lembrar que dentro do contexto de negócio o processo de seleção de variáveis conta com o apoio dos especialistas da área, o que ajuda bastante, porém as vezes por questões intrínsecas do negócio uma variável que é rejeitada por algum método de seleção acaba permanecendo no modelo.

fit.model = glm(SalePrice ~ MSSubClass + LotArea + OverallQual + 
    OverallCond + YearBuilt + YearRemodAdd + MasVnrArea + BsmtUnfSF + 
    TotalBsmtSF + X2ndFlrSF + GrLivArea + BsmtFullBath + BedroomAbvGr + 
    KitchenAbvGr + TotRmsAbvGrd + Fireplaces + GarageYrBlt + 
    GarageCars + GarageArea + WoodDeckSF + MSZoning + Neighborhood +  
    HouseStyle + Exterior1st + MasVnrType + ExterQual + BldgType + 
    Foundation + BsmtQual + BsmtExposure + BsmtFinType1 + HeatingQC + 
    CentralAir + KitchenQual + Functional + SaleCondition,
    data = model.train,family=Gamma(link=identity))

summary(fit.model)
## 
## Call:
## glm(formula = SalePrice ~ MSSubClass + LotArea + OverallQual + 
##     OverallCond + YearBuilt + YearRemodAdd + MasVnrArea + BsmtUnfSF + 
##     TotalBsmtSF + X2ndFlrSF + GrLivArea + BsmtFullBath + BedroomAbvGr + 
##     KitchenAbvGr + TotRmsAbvGrd + Fireplaces + GarageYrBlt + 
##     GarageCars + GarageArea + WoodDeckSF + MSZoning + Neighborhood + 
##     HouseStyle + Exterior1st + MasVnrType + ExterQual + BldgType + 
##     Foundation + BsmtQual + BsmtExposure + BsmtFinType1 + HeatingQC + 
##     CentralAir + KitchenQual + Functional + SaleCondition, family = Gamma(link = identity), 
##     data = model.train)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.95166  -0.05437   0.00005   0.05192   0.40735  
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -8.441e+05  1.166e+05  -7.238 8.15e-13 ***
## MSSubClass           -1.158e+02  5.175e+01  -2.238 0.025386 *  
## LotArea               1.066e+00  1.724e-01   6.186 8.48e-10 ***
## OverallQual           5.828e+03  7.163e+02   8.136 1.01e-15 ***
## OverallCond           4.953e+03  5.526e+02   8.963  < 2e-16 ***
## YearBuilt             2.999e+02  4.846e+01   6.187 8.41e-10 ***
## YearRemodAdd          6.302e+01  3.576e+01   1.762 0.078250 .  
## MasVnrArea            8.576e+00  5.319e+00   1.612 0.107147    
## BsmtUnfSF            -8.395e+00  2.394e+00  -3.507 0.000471 ***
## TotalBsmtSF           2.694e+01  3.961e+00   6.801 1.64e-11 ***
## X2ndFlrSF             1.753e+01  5.209e+00   3.365 0.000791 ***
## GrLivArea             4.510e+01  4.026e+00  11.202  < 2e-16 ***
## BsmtFullBath          2.043e+03  1.281e+03   1.595 0.111087    
## BedroomAbvGr         -3.979e+03  9.753e+02  -4.080 4.79e-05 ***
## KitchenAbvGr         -2.069e+04  3.885e+03  -5.327 1.19e-07 ***
## TotRmsAbvGrd          1.855e+03  7.252e+02   2.558 0.010640 *  
## Fireplaces            2.953e+03  9.755e+02   3.027 0.002520 ** 
## GarageYrBlt           7.349e+01  3.656e+01   2.010 0.044637 *  
## GarageCars            3.583e+03  1.578e+03   2.271 0.023344 *  
## GarageArea            1.119e+01  5.874e+00   1.905 0.057027 .  
## WoodDeckSF            1.118e+01  4.520e+00   2.474 0.013489 *  
## MSZoningFV            4.486e+04  8.240e+03   5.445 6.30e-08 ***
## MSZoningRH            3.569e+04  6.307e+03   5.659 1.90e-08 ***
## MSZoningRL            3.813e+04  4.721e+03   8.077 1.61e-15 ***
## MSZoningRM            3.301e+04  4.107e+03   8.038 2.18e-15 ***
## NeighborhoodBlueste   1.952e+03  1.211e+04   0.161 0.871949    
## NeighborhoodBrDale    1.163e+04  7.605e+03   1.529 0.126516    
## NeighborhoodBrkSide   6.684e+03  7.113e+03   0.940 0.347590    
## NeighborhoodClearCr  -6.309e+03  7.791e+03  -0.810 0.418212    
## NeighborhoodCollgCr  -5.532e+03  5.999e+03  -0.922 0.356623    
## NeighborhoodCrawfor   2.048e+04  6.862e+03   2.985 0.002893 ** 
## NeighborhoodEdwards  -1.122e+04  6.348e+03  -1.767 0.077519 .  
## NeighborhoodGilbert  -6.596e+03  6.425e+03  -1.027 0.304795    
## NeighborhoodIDOTRR    3.103e+03  7.630e+03   0.407 0.684285    
## NeighborhoodMeadowV  -4.235e+03  8.621e+03  -0.491 0.623340    
## NeighborhoodMitchel  -1.026e+04  6.615e+03  -1.551 0.121120    
## NeighborhoodNAmes    -9.638e+03  6.202e+03  -1.554 0.120452    
## NeighborhoodNoRidge   3.150e+04  8.088e+03   3.895 0.000104 ***
## NeighborhoodNPkVill   1.112e+04  7.918e+03   1.405 0.160305    
## NeighborhoodNridgHt   1.841e+04  6.706e+03   2.745 0.006148 ** 
## NeighborhoodNWAmes   -1.001e+04  6.490e+03  -1.543 0.123171    
## NeighborhoodOldTown  -4.863e+03  7.151e+03  -0.680 0.496662    
## NeighborhoodSawyer   -6.549e+03  6.418e+03  -1.020 0.307786    
## NeighborhoodSawyerW  -3.635e+03  6.381e+03  -0.570 0.568993    
## NeighborhoodSomerst   5.286e+03  8.398e+03   0.629 0.529185    
## NeighborhoodStoneBr   3.403e+04  8.063e+03   4.220 2.63e-05 ***
## NeighborhoodSWISU     5.289e+03  7.194e+03   0.735 0.462392    
## NeighborhoodTimber   -7.486e+03  7.166e+03  -1.045 0.296382    
## NeighborhoodVeenker   3.836e+03  9.242e+03   0.415 0.678159    
## HouseStyle1.5Unf      1.125e+04  5.059e+03   2.223 0.026400 *  
## HouseStyle1Story      8.747e+03  2.829e+03   3.092 0.002034 ** 
## HouseStyle2.5Unf      6.750e+03  5.962e+03   1.132 0.257769    
## HouseStyle2Story     -6.964e+02  2.345e+03  -0.297 0.766566    
## HouseStyleSFoyer      1.105e+04  4.368e+03   2.529 0.011563 *  
## HouseStyleSLvl        1.126e+04  3.737e+03   3.012 0.002646 ** 
## Exterior1stBrkComm   -4.424e+04  9.799e+03  -4.515 6.96e-06 ***
## Exterior1stBrkFace    3.063e+03  4.264e+03   0.718 0.472777    
## Exterior1stCemntBd   -4.334e+03  5.919e+03  -0.732 0.464180    
## Exterior1stHdBoard   -1.072e+04  3.661e+03  -2.928 0.003477 ** 
## Exterior1stMetalSd   -4.994e+03  3.448e+03  -1.448 0.147767    
## Exterior1stPlywood   -1.039e+04  4.110e+03  -2.528 0.011585 *  
## Exterior1stStucco    -3.908e+03  4.929e+03  -0.793 0.427954    
## Exterior1stVinylSd   -7.128e+03  3.542e+03  -2.013 0.044382 *  
## Exterior1stWd Sdng   -8.151e+03  3.516e+03  -2.318 0.020596 *  
## Exterior1stWdShing   -9.666e+03  4.692e+03  -2.060 0.039603 *  
## MasVnrTypeBrkFace     8.424e+03  4.145e+03   2.032 0.042330 *  
## MasVnrTypeNone        9.125e+03  4.195e+03   2.175 0.029802 *  
## MasVnrTypeStone       1.269e+04  4.768e+03   2.661 0.007887 ** 
## ExterQualFa          -3.218e+04  7.953e+03  -4.046 5.55e-05 ***
## ExterQualGd          -2.008e+04  6.111e+03  -3.286 0.001046 ** 
## ExterQualTA          -2.379e+04  6.284e+03  -3.785 0.000161 ***
## BldgType2fmCon        1.416e+04  7.529e+03   1.881 0.060277 .  
## BldgTypeDuplex       -6.343e+03  4.889e+03  -1.297 0.194725    
## BldgTypeTwnhs        -1.342e+04  6.609e+03  -2.030 0.042536 *  
## BldgTypeTwnhsE       -4.731e+03  5.954e+03  -0.795 0.427017    
## FoundationCBlock      6.254e+03  1.972e+03   3.171 0.001556 ** 
## FoundationPConc       6.202e+03  2.196e+03   2.825 0.004808 ** 
## FoundationStone       1.810e+04  6.244e+03   2.899 0.003808 ** 
## FoundationWood       -1.817e+04  1.108e+04  -1.640 0.101208    
## BsmtQualFa           -2.224e+04  4.616e+03  -4.818 1.64e-06 ***
## BsmtQualGd           -2.166e+04  3.519e+03  -6.156 1.02e-09 ***
## BsmtQualTA           -2.309e+04  3.890e+03  -5.937 3.79e-09 ***
## BsmtExposureGd        1.058e+04  2.595e+03   4.078 4.84e-05 ***
## BsmtExposureMn       -3.661e+03  2.417e+03  -1.514 0.130172    
## BsmtExposureNo       -3.650e+03  1.846e+03  -1.978 0.048197 *  
## BsmtFinType1BLQ       1.567e+02  1.796e+03   0.087 0.930457    
## BsmtFinType1GLQ       3.790e+03  1.887e+03   2.009 0.044779 *  
## BsmtFinType1LwQ      -5.350e+03  2.323e+03  -2.303 0.021476 *  
## BsmtFinType1Rec      -1.618e+03  1.888e+03  -0.857 0.391794    
## BsmtFinType1Unf      -3.108e+03  2.016e+03  -1.542 0.123407    
## HeatingQCFa           3.052e+03  3.008e+03   1.014 0.310573    
## HeatingQCGd          -1.430e+03  1.459e+03  -0.980 0.327152    
## HeatingQCPo          -2.351e+04  1.097e+04  -2.143 0.032317 *  
## HeatingQCTA          -2.675e+03  1.382e+03  -1.935 0.053188 .  
## CentralAirY           3.736e+03  2.141e+03   1.745 0.081258 .  
## KitchenQualFa        -2.383e+04  4.824e+03  -4.940 8.93e-07 ***
## KitchenQualGd        -2.266e+04  3.654e+03  -6.203 7.64e-10 ***
## KitchenQualTA        -2.244e+04  3.756e+03  -5.975 3.03e-09 ***
## FunctionalMaj2       -4.983e+03  7.777e+03  -0.641 0.521827    
## FunctionalMin1        5.310e+03  6.313e+03   0.841 0.400463    
## FunctionalMin2        2.512e+03  6.218e+03   0.404 0.686356    
## FunctionalMod         7.177e+02  8.199e+03   0.088 0.930256    
## FunctionalTyp         1.384e+04  5.589e+03   2.476 0.013408 *  
## SaleConditionAdjLand  2.315e+04  1.383e+04   1.674 0.094376 .  
## SaleConditionAlloca  -3.831e+03  6.013e+03  -0.637 0.524144    
## SaleConditionFamily   2.856e+03  3.692e+03   0.773 0.439462    
## SaleConditionNormal   3.334e+03  1.703e+03   1.958 0.050430 .  
## SaleConditionPartial  1.120e+04  3.161e+03   3.543 0.000411 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01073054)
## 
##     Null deviance: 189.616  on 1299  degrees of freedom
## Residual deviance:  13.611  on 1192  degrees of freedom
## AIC: 29343
## 
## Number of Fisher Scoring iterations: 7

3.1 Teste Qui-Quadrado

round(1-pchisq(13.611,1192, ncp = 0, lower.tail = TRUE, log.p = FALSE),3)
## [1] 1

O desvio do modelo foi de D(\(y\);\(\hat{\mu}\)) = 13.611, com 1192 graus de liberdade, que leva a P = 1,00 e indica um ajuste adequado.

4 Dados de teste

model.test <- fread('../outputs/model.test.csv', 
                     sep=",", 
                     showProgress = FALSE)[,-1] %>%
                     data.frame(stringsAsFactors = T) %>% 
                     select(Id,everything())  
tipo <- lapply(model.test,class)
model.test[,unlist(tipo) != 'integer'] <- data.frame(apply(model.test[,unlist(tipo)!='integer'],2,factor))
model.test

4.1 Predição

predicao <- predict(fit.model,model.test)

4.2 Plot Ecdf

Como não temos os valores da variável SalePrice no conjunto de dados de teste, podemos apenas verificar se a distribuição dos dados preditos esta de acordo com a distribuição dos dados reais da variável SalePrice do conjunto de dados de treino. A da distribuição impírica acumulada dos dados pode nos ajudar com isso.

latticeExtra::ecdfplot(~ fit.model$data$SalePrice + predicao, 
                         auto.key=list(space='bottom',col = c('red','blue')), 
                         col = c('red','blue'), 
                         lwd = c(2,3), 
                         xlab =" ",ylab = 'F(x)',
                         main = 'Distribuição Empírica Acumulada')

Observe que há concordância entre os valores preditos e reais da variável SalePrice. Uma outra forma de verificar essa condordância é plotar as distribuições de densidade e frequência, como fiz abaixo.

predicao <- data.frame(predicao)
p1 <- ggpubr::ggdensity(predicao, x = "predicao", 
              fill = "#0073C2FF", color = "black",
              add = "mean", rug = TRUE) +
              xlim(xlim = c(0,max(model.train$SalePrice)) ) +
              ylim(ylim = c(0,8.0*10^(-6))) +    
              labs(title = 'Distribuição da densidade das predições') +
              theme_dark()

p2 <- ggpubr::ggdensity(model.train, x = "SalePrice", 
              fill = "lightyellow", color = "black",
              add = "mean", rug = TRUE) +
              xlim(xlim = c(0,max(model.train$SalePrice)) ) +
              ylim(ylim = c(0,8.0*10^(-6) )) +  
              labs(title = 'Distribuição da densidade da variável SalePrice') +
              theme_dark()

gridExtra::grid.arrange(p1,p2,nrow = 1)

Embora as distribuições estejam muito parecidas, os gráficos acima nos permite visualizar que os valores preditos assume uma variabilidade um pouco maior a partir do valor médio da distribuição.

p3 <- ggplot(predicao, aes(x = predicao)) + 
      geom_histogram(bins = 100, 
                     color = "black", 
                     fill = "#0073C2FF") +
      xlim(xlim = c(0,max(model.train$SalePrice)) ) +
      ylim(ylim = c(0,100)) +
      labs(title = 'Distribuição de frequência das predições') +
      theme_dark()

p4 <- ggplot(model.train, aes(x = SalePrice)) + 
      geom_histogram(bins = 100, 
                     color = "black", 
                     fill = "lightyellow") +
      xlim(xlim = c(0,max(model.train$SalePrice)) ) +
      ylim(ylim = c(0,100)) +
      labs(title = 'Distribuição de frequência da variável SalePrice') +
      theme_dark() 


gridExtra::grid.arrange(p3,p4,nrow = 1)

Se no conjunto de teste a variável SalePrice mantiver a mesma distribuição do conjunto de treino meu modelo terá boas chances de performar bem.

plot(fit.model$data$SalePrice,
     fit.model$fitted.values,
     col = c('blue','red'), 
     pch = 20,
     type = 'p',
     xlab = 'SalePrice',
     ylab = 'Valores Ajustados',
     main = 'SalesPrice vs Valores Preditos')
abline(0,1)
legend('topleft',legend = c('SalePrice','Valores ajustados'),col = c('blue','red'), lty = c(2,2) )

Este plot mostra que nosso modelo é eficaz para valores até \(4e+05\), isso ocorre por existirem poucas casas com valores muito alto no conjunto de dados, pois a variabilidade tem convergência assintótica, ou seja, \(\sigma_{_{n \implies \infty}} \implies 0\). Uma alternativa visando melhorar a performance do modelo para valores acima de \(4e+05\) é criar variáveis que ajudam a discriminar a variável resposta neste segmento. Através de uma clusterização acho que isso seria possível.

4.3 Resíduos

plot(fit.model$residuals, pch = 20, type = 'p', main = 'Resíduos')

Os resíduos estão bem distribuídos em torno de zero e isso é muito bom! Considerando que o ajuste está bom, é interessante notar que os resíduos que estão muito acima de zero podem ser casos em que o valor do imóvel esta fora do praticado pelo mercado e isto implica que dificilmente um imóvel nesta situação seria facilmente vendido fazendo com que a equipe de vendas perdesse tempo e energia. Já os imóveis com preços abaixo de zero seriam facilmente vendidos porém o lucro da venda seria baixo.

5 Rsquare

rss <- sum((fit.model$residuals)^2)  ## soma dos quadrados dos resíduos 
tss <- sum((fit.model$data$SalePrice - mean(fit.model$data$SalePrice))^2)  ## soma total dos quadrados
(r.square <- 1-rss/tss)
## [1] 0.8942537

Este é um bom resultado para uma primeira iteração, o modelo consegue explicar cerca de 89% da variabilidade dos dados.

6 RMSE

sqrt(mean((fit.model$residuals)^2))
## [1] 25378.23

Com esse RMSE somente o kaggle me dirá se está bom ou não.

7 Gráficos de Dispersão

layout(matrix(c(1,1,2,2,3,3), ncol = 2, nrow = 3, byrow = T))
plot(fit.model$data$SalePrice,
     col = 'black', 
     pch = 20,
     type = 'p',
     ylim = c(0,7*10^5),
     xlab = 'SalePrice',
     ylab = '',
     main = 'Valores reais de SalesPrice')
plot(fit.model$fitted.values,
     col = 'blue', 
     pch = 20,
     type = 'p',
     ylim = c(0,7*10^5),     
     xlab = 'SalePrice',
     ylab = '',
     main = 'Valores Ajustados de SalesPrice')
plot(predicao$predicao,
     col = 'red', 
     pch = 20,
     type = 'p',
     ylim = c(0,7*10^5),     
     xlab = 'SalePrice',
     ylab = '',
     main = 'Valores Preditos com df.test para SalesPrice')

8 QQ-Plot

O QQ-plot pode nos ajudar a ter uma idéia de como será o resultados da predição, aparemtemente será bom.

qqplot(sort(fit.model$data$SalePrice)[1:1318],
       sort(predicao$predicao),
       main = 'QQ-plot',
       ylab = 'Predição',
       xlab = 'SalePrice',
       col = c('blue','red') )
abline(0,1)

9 Exportando para submissão no Kaggle

submit.glm <- data.frame(Id = model.test$Id,SalePrice = predicao$predicao, stringsAsFactors = F)
write.csv(submit.glm, file="../outputs/saleprice_glm.csv",row.names = FALSE,quote=FALSE)

10 Conclusão

Neste primeiro momento os resultados são razoáveis.

Sérgio Carvalho

20 março, 2019